Phase field modeling of electrochemistry I: Equilibrium 
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A diffuse interface (phase field) model for an electrochemical system is developed. We describe 
the minimal set of components needed to model an electrochemical interface and present a varia- 
tional derivation of the governing equations. With a simple set of assumptions: mass and volume 
constraints, Poisson's equation, ideal solution thermodynamics in the bulk, and a simple description 
of the competing energies in the interface, the model captures the charge separation associated with 
the equilibrium double layer at the electrochemical interface. The decay of the electrostatic potential 
in the electrolyte agrees with the classical Gouy-Chapman and Debye-Hiickel theories. We calculate 
the surface free energy, surface charge, and differential capacitance as functions of potential and find 
qualitative agreement between the model and existing theories and experiments. In particular, the 
differential capacitance curves exhibit complex shapes with multiple extrema, as exhibited in many 
electrochemical systems. 



I. INTRODUCTION 

We develop a phase field model of an electrochemi- 
cal system. The method employs a phase field variable, 
which is a function of position and time, to describe 
whether the material is one phase or another (i.e., the 
electrode or electrolyte). The behavior of this variable is 
governed by a partial differential equation (PDE) that is 
coupled to the relevant transport equations for the ma- 
terial. The interface between the phases is described by 
smooth but highly localized changes of this variable. This 
approach avoids the mathematically difficult problem of 
applying boundary conditions at an interface whose lo- 
cation is part of the unknown solution. The phase field 
method is powerful because it easily treats complex in- 
terface shapes and topology changes. One long range 
goal of the approach is to treat the complex geometry, 
including void formation, that occurs during plating in 
vias and trenches for on-chip metallization [1]. 

Early models of the electrochemical interface, devel- 
oped by Gouy, Chapman, and Stern [2, 3, 4, 5] focused 
on the distribution of charges in the electrolyte. These 
models, which assume that the charges have a Boltzmann 
distribution and are subject to Poisson's equation, arc 
summarized in Appendix D. More recently, density func- 
tional models have been applied to the equilibrium elec- 
trochemical interface [6, 7]. These atomic scale models 
describe the electrolyte with distribution functions which 
have maxima at the positions of the atoms and take the 
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electrode to be a hard, idealized surface. The equilib- 
rium distribution of electrons and ions are computed in 
these models and relationships between potential, charge, 
surface free energy, and capacitance are obtained. No ki- 
netic modelling is performed in these papers. Phase field 
models can be viewed as a mean field approximation of 
atomic scale density functional theories [8, 9, 10], and the 
two methods often make similar predictions. 

The phase field method has been used widely for solid- 
ification [8, 9]. The present approach is motivated by the 
mathematical analogy between the governing equations 
of solidification dynamics and electroplating dynamics. 
For example, the solid-melt interface is analogous to 
the electrode-electrolyte interface. The various ovcrpo- 
tentials of electrochemistry have analogies with the su- 
percoolings of alloy solidification: diffusional (constitu- 
tional), curvature and interface attachment. Dendrites 
can form during solidification and during electroplating. 
It is not surprising, however, that we find significant dif- 
ferences between the two systems. The crucial presence 
of charged species in electrochemistry leads to rich in- 
teractions between concentration, electrostatic potential, 
and phase stability. 

We first pick a minimal set of components required 
to describe the possible composition variations from the 
electrode to the electrolyte through the electrochemical 
interface. Next we define concentration and mole fraction 
variables. Then a variational principle is used to establish 
a set of PDE's that govern equilibrium interfaces. In a 
second paper [11], we explore dynamic solutions to the 
phase field equations. 

Given a minimal set of assumptions, the model predicts 
the charge separation associated with the equilibrium 
double layer at the electrochemical interface. Changes 
in surface potential induce changes in surface free en- 
ergy (electrocapillary curves), surface charge, and differ- 
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cntial capacitance. The decay length of charge in the 
electrolyte as a function of electrolyte concentration is 
consistent with the Debye-Hiickel theory. The parame- 
ters of the phase field model are related to the physical 
parameters of traditional electrochemistry. The results 
are compared to the classic Gouy-Chapman-Stern model 
of the electrified interface as well as to experimental ca- 
pacitance curves [12, 13]. 

II. MODEL FORMULATION 



the sum over all substitutional components. The mole 

per volume concentrations are defined as Cj = Xj/V m , 
such that 

n n 

B C i = ® = L ( 4 ) 

3=1 j=2 

Picking one component j = n with V n = V s ^ 0, we can 
always express its concentration in terms of the others as 

c„ = i/y s -E"=2^. 



A. Choice of Components, Phases, and Molar 
Volumes 

To treat an electrochemical system, we consider a set of 
components consisting of cations, anions, and electrons. 
We will refer to the metallic conducting electrode phase 
as a and the ionic conducting electrolyte phase as (3. 
In this diffuse interface model, the concentrations of the 
components will vary smoothly through the interfacial re- 
gion. The electrode will be considered a phase of primar- 
ily two components: electrons e - (component #1) and 
cations M +m (component #2). The electrolyte will be 
considered a phase of primarily three components; more 
noble cations M +m (component #2), less noble cations 
N + ™ (component #3), and anions A~ a (component #4). 
A model aqueous electrolyte can be considered as a spe- 
cial case by setting n = 0. The primary charge transfer 
reaction for this system is 

M +Tn (/?) + m e" (a) =± M (a) . (1) 

All ions are treated as substitutional species with iden- 
tical partial molar volumes V s . Electrons are treated as 
an interstitial species with zero partial molar volume. 
The interstitial nature is necessary to recover ohmic be- 
havior in the electrode phase. In reality, the partial molar 
volumes of the substitutional components should depend 
at least on the species and phase, but this variation leads 
to deformation and flow, which are not the focus of this 
work. We employ mole fractions Xj of each component 
j that have the conventional definition, such that 

n 
3=1 

where n = 4 for the four-component system we consider 
in this paper. The molar volume varies because the mole 
fractions of the components vary from one phase, through 
the interface, and into the other phase. At each point it 
is given by 

n n 

V m = Y,VjX j = V s Y J X h (3) 

3=1 3=2 

where Vj is the partial molar volume of each component 
j, Yl^—i is the sum over all components, and E" =2 is 



B. Equilibrium 

We propose a Helmholtz free energy for an isothermal 
system of charged components 

F C n , ft 

= f v (fv{Z,C 1 ...C n ) + ^\S7tf + ±p<f?\dV (5) 

integrated over the volume V, where fv is the Helmholtz 
free energy per unit volume, £ is the phase field variable, 
4> is the electrostatic potential, is the gradient energy 
coefficient for the phase field, 

pz-F^zM (6) 

3=1 

is the charge density, Zj is the valence of component j 
(equiv/mol or charge units per ion), and T is Faraday's 
constant. The first term in the integral of Eq. (5) rep- 
resents the energy density of a system without gradients 
and with no charge interactions; the second represents 
the gradient energy without electrostatic effects; and the 
third represents the electrostatic energy. 

There are three constraints on the field variables of this 
system. The total number of moles rij of each species j 
must be conserved over the volume V 

rij = [ Cj dV = VCj, j = l...n (7) 
Jv 

where Cj is the average concentration of species j. In 
addition, Eq. (4) and Poisson's equation 

V • [e(£)V0] + p = (8) 

must be satisfied at every point in the system. e(£) is the 
electrical permittivity, whose value we take to explicitly 
depend on the phase. Since the phase field is coupled 
to the other variables, there is also an implicit depen- 
dence on the concentrations of different species in the 
electrolyte. 

We note that by invoking Poisson's equation (8), the 
electrostatic energy term in Eq. (5), given by ^p4>, could 
cquivalently be represented as • E, one half the scalar 
product of the displacement and the electric field of elec- 
tromagnetic theory, or as ||V0| 2 . In this third form, 
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we see that the electrostatic energy density is completely 
analogous to the phase field gradient energy density. 

We perform a variational analysis on the resulting La- 
grangian 

C = F - jf Ay(x) \J2 VjCj - lj dV 

n 



i=i 



A (x){V- [e(0V4>]+p}dV (9) 



where we have introduced the Lagrange multipliers 
Ay(x), Xj, and A^x) for the constraints (4), (7), and 
(8) (note that the requirement that Eqs. (4) and (8) be 
satisfied everywhere in the system means that Ay and A^ 
must be fields). At equilibrium, each of the variations of 
C must be independently zero. 

= = §FT + \^ Z A - h ~ Fzj\,p - Ay Vj (10a) 



3 = 1 



^ = O=^-- K? V 2 e + e'(OVA V0 
g = 0=ip-V-[e(OVA,]. 



(10b) 
(10c) 



Using Poisson's equation to eliminate p in Eq. (10c), we 
determine that = —<f)/2. Here we have assumed that 
the natural boundary conditions = d£,/dn = d(j>/dn = 
dX^/dn hold on the boundary of V. Substituting this 
result into Eq. (10a), 



+ Tz$- \v(x)Vj. j = l. 



(11) 



Since V e - = and Vj = V s for all substitutional species, 
we can eliminate Ay(x). 

We thus can summarize the equilibrium governing 
equations as 



Xjn 



A, - ^A r 



df 



v 



= 



dCi 
constant 
dfv 



J-z 



V, 



dfv 



= V • [e(OV0] 4 



dC n 
j = 1 ...» - 1 

2 C 



■ J~ z n 



(12a) 
(12b) 
(12c) 



and the classical electrochemical potentials 
dfv 



da 1 - 



j = l...n (14) 



such that Eq. (12a) becomes 
Vj _ 

[ij — ryi^n = constant. 
In Appendix A we show that, in one dimension, 



(15) 



fh=Pr + Vj^-^), J = l...n (16) 

where £l°° is the electrochemical potential far from the in- 
terface. First, we note that far from the interface, where 
V£ = V</> = 0, the value of the electrochemical poten- 
tials for each component j are identical in each phase, in 
agreement with Gibbsean thermodynamics [14] 



p5- 



3 = 1 



(17) 



Second, the electrochemical potential of e _ is uniform 
throughout the system. For the substitutional compo- 
nents, the electrochemical potential varies through the 
interface, even though the values in the bulk phases are 
equal. Because = far from the interface, in the ab- 
sence of external charges, Eq. (12c) requires that charge 
is zero in the bulk phases, such that 



Y Z 3 X j - Z 3 X J 



(18) 



III. INTERFACIAL PROPERTIES 

From electrocapillary theory [2] , we know that surface 
free energy 7, excess charge on the electrode a a . and 
differential capacitance Cd are related by 



( dl \ 



( d 2 7 



(19) 
(20) 



where = <f) a — <fiP is the applied potential difference 
across the interface. 1 The potential in the electrolyte 
far from the interface is and that in the electrode is 
4> a . In a perfect conductor is uniform throughout 
the phase. The variation at constant chemical potential 



It is convenient to identify the classical chemical po- 
tentials 



df 



v 



3 = 1 . . . n 



(13) 



1 In classical electrochemistry, the actual difference in electrostatic 
potential across the interface is called the "inner" or "Galvani" 
potential difference. 
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(constant concentration) is difficult to perform with phys- 
ical electrode/electrolyte systems. In general, this vari- 
ation can only be performed experimentally for an inert 
system such as a mercury electrode against an aqueous 
electrolyte. We seek to determine whether these relation- 
ships hold for our phase field model. 

In Appendix B, we derive the expression for the surface 
free energy 



(21) 



The first term in the integrand represents the contribu- 
tions of the non-electrified interface. The second term is 
the contribution of electrostatics, such that the presence 
of electric fields always reduces the surface free energy 
from its charge-free value. In Appendix C we obtain the 
relationship 



\dA(f>° 



if <7 Q is defined by 



p{C)pdx, 



(22) 



(23) 



which is a completely reasonable definition of the surface 
charge on the electrode. p(£) is an interpolation function 
that will be described in Section IV A; in short, = 1 
in the electrode and = in the electrolyte. The no- 
tation A(f>° differs from that in Eqs. (19) and (20). A(j)° 
refers to a materials property of the electrode-electrolyte 
system, as we will discuss in Section IV B. Because we 
consider a non-inert electrode, we cannot vary the ap- 
plied potential without affecting the concentration in the 
electrolyte. The variation considered in Eq. (22) is ac- 
tually a variation with respect to a changing material 
property, rather than an applied potential. The differen- 
tial capacitance is then defined to be 



a 



dc 



\dA4> c 



(24) 



We note that a curved interface exhibits an additional re- 
lationship between surface free energy, surface concentra- 
tion, and the interfacial potential drop, consistent with 
the Gibbs-Thomson effect [15]. 



A. Choice of Form of the Thermodynamic Function 

For simplicity, we assume that the chemical part of the 
Helmholtz free energy per unit volume is described by an 
interpolation of two ideal solutions of the components for 
the electrode and electrolyte. 

fv (£j C\... C n ) = —fm (£, X\ . . . X n ) 

n 

3=1 

+ RT\nC j V m + W j g{^)} (25) 

where and fi°^ are the chemical potentials of pure 
component j in the electrode (metal) phase and the elec- 
trolyte phase respectively, R is the molar gas constant, 
and T is the temperature. Following many phase field 
models of solidification [16], we use an interpolating func- 
tionp (£) = £ 3 (6£ 2 - 15£ + 10) to bridge between the de- 
scriptions of the two bulk phases and a double- well func- 
tion g (£) = ^ 2 (1 — with a barrier height of Wj for 
each component j to establish the metal/electrolyte in- 
terface. The barrier heights Wj penalize interfaces which 
are too broad and the gradient energy coefficient pe- 
nalizes interfaces which arc too narrow. The polynomials 
are chosen to have the properties that p(Q) = 0, p(i) = 1, 
2/(0) = 2/(1) = 0, and g'(0) = g'(l) = 0. Other functions 
could be used. 

For use in Eq. (12b), this free energy leads to 



dfy 



p' (0 ]T CjAfij + g' (0 ]T CjWj, (26) 



where Atf = - nf ■ 

The quantity dfv/dCj is also computed to give the 
electrochemical potentials, 

flj = nf + AfijP (0 + RT In CjV m + ZjT<f> + W j9 (0 . 

j = 1. . .n (27) 

We note that the p,j depend on all the Cj through the 
molar volume V m . 



IV. THERMODYNAMIC FUNCTIONS AND 
MATERIAL PARAMETERS 

To make the results of the phase field model more 
concrete and to permit numerical calculations, we must 
choose a particular form of the thermodynamic function 
and the materials parameters, at which point the model 
will be fully specified. 



B. Standard Chemical Potentials 

We require values for the parameters A/j,° for each of 
the n species in our model in order to perform numer- 
ical calculations with the ideal solution thermodynamic 
model employed here. Given these numbers, we will know 
how the bulk electrode and electrolyte concentrations 
vary with the potential difference across the interface. 
From the equality of the bulk electrochemical potentials 
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Eq. (17) and the ideal solution form of the electrochem- 
ical potential Eq. (27), 



Acj) 



j = 1 . . . n 



(28) 



where the mole fractions are constrained by charge neu- 
trality of the bulk phases Eq. (18) and the ordinary def- 
inition of the mole fractions Eq. (2). 

One procedure to obtain information about the A/x? 
is to specify one set of concentrations that are at equi- 
librium, X"° and Xj° . With such a set we can only de- 
termine three linear combinations of A/j,°. Given these 
three linear combinations, we cannot determine the po- 
tential difference across the interface, but we can calcu- 
late how the bulk electrode and electrolyte concentrations 
vary with changes in the potential difference across the 
interface. In other words, the potential difference can 
only be described with respect to a reference electrode. 
Although knowledge of the potential difference between 
the phases is not necessary to describe the bulk equilib- 
rium, it is necessary if one is interested in modeling the 
charge distribution between the electrode and electrolyte. 
We will designate A(j)° as the potential difference across 
the interface for the mole fractions X"° and X? , such 
that 



A M ° 



X 1 

AT In J 



0° 



3 = 1 



(29) 



Thus all four values of A/i° can be computed. The quan- 
tity A<f>° is a material property that is fixed for a given 
choice of electrode and electrolyte system. 

In traditional electrochemistry, one takes the reference 
mole fractions of the clcctroinactive species to be zero in 
one phase or the other, e.g., = X^l = 0. One would 
then only equate the electrochemical potentials between 
the bulk phases of the electroactive species M +m and 
N + ". In this case it can be shown that 



£- 



N+" 



mj 7 



(30) 



where the standard potentials are ^n+™ are 

tained from a table of electromotive series. Eq. (30) 
would be adequate to describe the concentration vari- 
ations between the bulk phases. 

For numerical purposes, we must assume small, but 
non-zero values for the reference concentrations of the 
electroinactive species, equivalent to assuming large (pos- 
itive or negative), but finite, values for their standard 
potentials. In reality these concentrations are not zero, 
as might be indicated by an electrolyte with electronic 
conductivity or an electrode with some anion solubility 
[2, 17, 18]. 

In the remainder of the paper, we perform a detailed 
equilibrium analysis for an electrolyte where the com- 
ponent N has no charge, like an aqueous system (N = 



TABLE I: Numerical values of the potential-independent por- 
tion of the chemical potential differences A/i° . 



in(xf /x;°) 


e 


-13.41 


M+ 2 


-2.919 


A" 2 


9.798 


N 


13.78 



electrode 
M 




MA 



electrolyte 



FIG. 1: Potential-composition phase diagram for the param- 
eters in Table I, illustrating the bulk equilibrium between a 
M electrode and a N electrolyte with dissolved MA salt. Tie- 
lines denote different values of the quantity {Acj) — A<j>°). The 
inset shows the position of this charge neutral phase diagram 
within the quaternary domain of the charged species. 



H2O) where dissociation of H2O is neglected. 2 The lower 
density of charged species in an unsupported "aqueous" 
electrolyte allows us to resolve the equilibrium interface 
more accurately than in a system where all species can 
carry charge. Our paper on the dynamic behavior of 
the electrochemical phase field model [11] will treat the 
case of N with charge. The bulk standard state mole 
fractions are chosen for an electrode of M metal and 
an electrolyte of a solution containing 1 mol/L of M +2 
and A -2 in a N solvent. We take the partial molar vol- 
ume of the "substitutional" components (M +2 , A -2 , and 
N) as that for pure water; i.e., V s = 0.018 L/mol or 
1/V S = 55.6 mol/L. The voltage-independent portion 
of the chemical potential differences are give in Table I. 



2 The dissociation of H2O could be handled by the addition of 
another component 
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The values for the electroinactive species are chosen to 
limit the corresponding standard state mole fractions to 
x 0o = x „o 2 = x «o = 1Q -6 

Bulk charge neutrality can be invoked to transform the 
four mole fraction variables X D - , X M +2, X A -2, and X-^ 
into mole fractions of three charge-neutral compounds, 
1m = |^ e -! ^ma = 2X A -2, and An- We plot the 
equilibrium phase diagram in terms of these transformed 
mole fraction coordinates in Figure 1. Equilibrium states 
exist only between -0.0898 V < A0 < +0.0427 V. It can 
be seen that the electrode remains essentially pure M 
(A^ +2 = 1/3, X"_ = 2/3) over the entire voltage range. 
At the lower A(f> limit, the electrolyte is essentially pure 
N (Xg = 1). At the upper A(j> limit, the electrolyte is 
pure MA (A^ +2 = 1/2, A^_ 2 = 1/2). 



V. NUMERICAL METHODS 

Numerical solutions to the governing equations 
Eq. (12) were obtained by both a relaxation method and 
a pseudo-spectral technique on a one-dimensional domain 
of length L. The relaxation method had the advantages 
that it was simple to code and would eventually converge 
to a solution even from a step-function initial condition. 
The main disadvantage of the technique is that it is very 
slow. Run times of several hours to several days were 
needed to reach convergence on a 1.6 GHz AMD Athlon 
running Dcbian GNU/Linux v3.0 with a 2.4 kernel and 
using the Portland Group pgee compiler. 3 In contrast, 
the pseudo-spectral technique can produce solutions in a 
few minutes, but only from a very good initial guess. 



Relaxation 



C. Other Parameters 

A simplification, which should be eliminated in future 
work, is employed for the barrier height Wj . We set equal 
values for the substitutional species Wj^2...n = W, and 
the value for the electrons, W e - = 0. This makes the 
last term in Eq. (26) independent of Cj and eliminates 
the Wj dependence in Eq. (12a). 

In phase field models of solidification, where electro- 
static effects arc not included, the surface free energy 7^ 
and the interfacial thickness 6% (of the £-field) are given 
by [19] 



7? = 

k = 



U 6 w 

2W ' 



(31) 
(32) 



We include the subscript £ with the realization that 
there will be an electrostatic contribution to the sur- 
face free energy (see Eq. (21)) and an independent elec- 
trostatic length scale. For numerical calculations, we 
choose the approximate values of 7^ = 0.2 J/m 2 and 
S s = 3 X 10" 11 m, which give W = 3.6 x 10 5 J/mol and 
= 3.6 x 10" 11 J/m. In reality, we would not expect 
a metallic electrode to have the same permittivity as an 
aqueous electrolyte. Moreover, the permittivity of the 
electrolyte is known to be lower near the interface than in 
the bulk as a result of the polarization in the double layer 
[4, 6]. While the variation in the permittivity undoubt- 
edly affects the structure of the interface, our goal in this 
paper is to show the richness obtained from a phase field 
model with even the simplest assumptions. We defer ex- 
amination of phase and concentration dependence of the 
permittivity to future work, and take e (£) = 78.49eo, 
where eo is the permittivity of free space. This is the 
value typically cited for an aqueous electrolyte [2]. 



Relaxation solutions [20] to Eq. (12) were obtained 
by casting the equilibrium ordinary differential equations 
(12a) and (12b) as the time dependent partial differential 
equations 



at 



and 

at 



5C 
5C~ 



I'j 



Yi 



(33a) 



3 = 1, 



-Mi 



SC 



-Me 



dfv 
OH 



(33b) 

where t is time, Mj is the mobility of component j, and 
is the mobility of the phase field. We defer discussion 
of the mobilities and Mj to our paper on kinetics [11] 
as their values are not important to the present analysis 
of electrochemical equilibrium. Equations (33) are the 
simplest expressions that guarantee a decrease in total 
free energy with time. Poisson's equation (12c) must still 
be satisfied everywhere. Equations (12c) and (33) were 
solved with explicit finite differences. Spatial derivatives 
were taken to second order on a uniform mesh. Solu- 
tions were integrated to equilibrium with an adaptive, 
fifth-order Rungc-Kutta time stepper [20]. We defined 
equilibrium as the point when Eqs. (12a) and (A6) were 
satisfied to within 0.1%. 

Simulations were started with an abrupt interface be- 
tween the bulk electrode and electrolyte phases, such that 



3 Certain commercial products are identified in this paper in order 
to adequately specify procedures being described. In no case 
does such identification imply recommendation or endorsement 
by the National Institute of Standards and Technology, nor does 
it imply the material identified is necessarily the best for the 
purpose. 
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TABLE II: Boundary Conditions 



electrode (x = 0) 


electrolyte (x — L) 


n ■ VC = 


n ■ V£ = 


<j> = 


n- = 


Cj specified 


Cj specified 



= 1 and £P = 0. After choosing a value for Cfl + , the 
remaining bulk Cj were determined from Figure 1. Be- 
cause p = for the bulk concentrations, Eq. (12c) gives 
V0 = throughout the domain for the initial data. The 
boundary conditions arc listed in Table II. 



(b) 





1.0- 
0.8 - 




■ 1 ■ 




0.6 - 
0.4- 
0.2 - 
0.0- 


electrode \ 


electrolyte 


(a) 










4- 






RT 






Adff/RT = 


0- 




A0°r/RT = 3.89 



B. Adaptive pseudospectral discretization 

In order to increase the numerical resolution of the 
intcrfacial region, we have also employed an adaptive so- 
lution technique based on a spectral approximation [21] 
to the governing equations (12). To reduce the number 
of unknowns in the system, we eliminate the solute vari- 
ables by solving the governing equations (12a), together 
with the constraint equation (4). Given values of £ and 
<j) at a point, this provides an algebraic form for the so- 
lute fields Cj at this point. The remaining equations 
(12b) and (12c) are discretized using a pseudospectral 
formulation of a spectral element representation of the 
second derivative [22]. It is convenient to fix the inter- 
face location by specifying = 1/2 at a given grid 
point xi. The discretization procedure then provides a 
set of nonlinear equations with an equal number of un- 
knowns. The nonlinear equations are solved using the 
quasi-Newton software package SNSQ [23]. Starting es- 
timates for the solution procedure are generally obtained 
by continuation from the finite difference procedure de- 
scribed above, or from previous pseudospectral solutions. 

An adaptive procedure is obtained by bisecting the el- 
ements for which an error estimate indicates that ad- 
ditional refinement is necessary. The error estimate is 
based on the rate of decay of a Chebyshev expansion of 
the solution components; a simple criterion is based on 
requiring the magnitude of the last two coefficients of the 
charge density in each panel to lie below a given thresh- 
old. If a refinement is necessary, the element is bisected 
and the previous solution is interpolated to the nodes of 
the new panels. The nonlinear equations are then solved 
on the new nodes, and the procedure is repeated until 
each Chebyshev expansion has rapid decay, indicating 
that the solution is well-resolved on each panel. Since 
the previous solution provides a good starting guess, the 
successive solutions converge quickly, and the overall run 
time is a small multiple of the time required for a single 
solution step on a grid with equal pseudospectral ele- 
ments. 



pVs 
T 



(c) 




x/L 



FIG. 2: Profiles through the interface of (a) the phase field 
variable, (b) the normalized electrostatic potential, and (c) 
the normalized charge distribution for A0° = -0.1 V, V, 
and +0.1 V. g(£) is mapped onto the background in gray to 
indicate the location of the phase field interface. The ticks 
marks at the top of Figure (a) indicate the positions of the 
mesh points for the two different solution methods. 



VI. NUMERICAL RESULTS 

Figure 2 shows plots of phase field, voltage and charge 
across the interface for A</>° = (-0.1, 0, +0.1) V obtained 
using the relaxation method in ID over a domain 3.2 nm 
long and containing 1200 mesh points. The tick marks at 
the top of Figure 2(a) indicate the positions of the mesh 
points used in the two different solution methods. The 
fields calculated by the relaxation and the pseudospec- 
tral methods are indistinguishable on the scale of these 
graphs, so we apply a linear correlation function to com- 
pare the calculation methods. The two methods have a 
linear correlation of 0.9992 for the most sensitive field 
p; the other fields have a linear correlation of 0.9999 or 
better. The difference between the two methods for the 
most sensitive field is thus of the same order as the cri- 
terion for stopping the relaxation calculations (< 0.1% 
error in Eqs. (12a) and (A6)); all other fields are much 
closer to convergence. 

Figure 2 focuses on the interface region of this compu- 
tational domain. The bulk concentration of M + and A" 
in the electrolyte is 0.25 mol/L. The variation of £ be- 
tween the electrode on the left and the electrolyte on the 
right for the three cases is virtually identical. A fit of all 
three curves to = {1 — tanh[(x — Xo)/26^]}/2 gives 
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FIG. 3: Normalized concentration profiles through the inter- 
face for different values of A(/}° . g(£) is mapped onto the 
background in gray to indicate the location of the phase field 
interface. 

5 e = (2.480 ± 0.009) x 10" 11 m, which compares well 
with the value we assumed in Section IV C. £ changes 
from 0.9 to 0.1 over a distance of approximately 0.1 nm 
or 4i$£. This represents the thickness of the electrode- 
electrolyte interface. The voltage <fi changes smoothly 
between a value of zero in the electrode (as assumed) to 
an asymptotic value in the electrolyte far from the inter- 
face equal to (RT/2F) ln[(l mol/L)/(0.25 mol/L)] - A</>°. 
The charge p, while zero far from the interfacial region, 
exhibits a distinct charge separation within the interfa- 
cial region. The charge distribution is quite different for 
the three cases. For A<fi° = -0.1 V and 0.0 V, a negative 
charge is present to the left. For A<j>° = +0.1 V, the 
negative charge is to the right. Figure 3 shows the varia- 
tion in the concentrations from the electrode to the elec- 
trolyte. The values at the ends of the full computational 
domain correspond to a bulk M electrode, a 0.25 mol/L 
MA in N electrolyte, with impurities as allowed by Ta- 
ble I. The abrupt change in concentrations through the 
distance where £ is changing is followed by a more grad- 
ual change in the electrolyte. The gradual concentration 
decay length in the electrolyte is the same as that of 
the voltage. One could define the surface excess as the 
difference between the actual concentration and some in- 
terpolation between the bulk values and see that there 
is an adsorption of the different species at the interface 
which depends on the value of A<fi°. 

In Appendix D we summarize the Gouy-Chapman- 
Stern model of the double layer. That treatment predicts 
an exponential decay of the potential in the electrolyte 
away from the electrode, with a decay length of <5? c . 




0.4 0.6 0.8 1.0 



x/L 

FIG. 4: Exponential fits (heavy dashed lines) to potential 
curves of Figure 2(b) (light solid lines). g(£) is mapped onto 
the background in gray to indicate the location of the phase 
field interface. 



Figure 4 shows a fit of the <j> vs. distance plots from Fig- 
ure 2(b) to (j) = 4>oo + (0o — 4>oo) exp(— x/S^). The fit 
is excellent. The decay length 8<j, of </> to its asymptotic 
value is very close to the predicted value of S^ c . This 
length is over ten times larger than 5^ and approximately 
three times the apparent interface thickness. 

Figure 5(a) shows the surface free energy (from 
Eq. (21)) vs. A(j)° and Figure 5(b) shows a plot of o a 
vs. A(jf , both obtained by the two numerical meth- 
ods. The point of zero charge (PZC), defined by a a — 
a' 3 = 0, occurs at A0° = +0.005 V for 1 mol/L and 
Acj) = +0.035 V for 0.25 mol/L. At the point of zero 
charge, p is not constant, nor is the electrostatic poten- 
tial, but rather, the integrated charge is zero in each 
phase and there is some potential step between them. 
We note that Grahame described exactly this condition 
in his seminal paper on the electrochemical double-layer 
[.]. The presence of dipoles at the interface guarantees 
that the potential will not be uniform. The surface charge 
curve shows a slight deviation from linearity away from 
the point of zero charge. This dome shaped curve in 
Figure 5(a) has a maximum surface free energy of ap- 
proximately 0.225 J/m 2 at a value of A0° = +0.005 V, 
the point of zero charge. This maximum surface free 
energy value is very close to 7^ , used to establish numer- 
ical values for W and k^. Figures 5(a) and 5(b) obey 
Eq. (19) very closely. The negative surface free energies 
obtained for large positive values of A<fi° indicate that 
a planar interface will become unstable to perturbations 
which increase surface area. Such perturbations are not 
possible given the symmetry constraints of our ID solu- 
tions, but attention will need to be paid to this when 
higher dimensional calculations are performed. 

From Eq. (20), the differential capacitance is obtained 
as the derivative of Figure 5(b) with respect to A0°. 
The relaxation method used to produce the open square 
points in Figure 5(b) is not fast to enough to allow cal- 
culating a numeric derivative of sufficient resolution. We 
thus use the results of the spectral method, which can 
compute with a much greater resolution and over a wider 
range of A(f>° , to calculate Figure 5(c). 

Our calculated differential capacitance curve, replot- 
ted in Figure 6(a), does not resemble the hyperbolic co- 
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FIG. 5: (a) Normalized surface free energy, (b) normalized 
surface charge, and (c) and normalized differential capaci- 
tance as functions of normalized A<j!> . Open symbols are cal- 
culated by the relaxation method of Section VA. Lines are 
calculated by the pseudospectral method of Section VB. 



sine predicted by the Gouy-Chapman theory {x-2 = in 
Eq. (D5)), shown in Figure 6(b). Neither does it resemble 
the truncated hyperbolic cosine predicted by the Gouy- 
Chapman-Stern theory (x2 ^ in Eq. (D5)), shown in 
Figure 6(c) (we take tjx-i — 5 F/m 2 for illustration only). 
On the other hand, it does bear a striking resemblance 
to experimental differential capacitance curves [2, 12, 13], 
such as Valette and Hamelin's measurements of Ag elec- 
trodes in NaF aqueous solutions, shown in Figure 6(d). 
The density functional calculations of Tang, Scriven, and 
Davis [7] also exhibit differential capacitance curves with 
multiple inflection points. 



VII. DISCUSSION AND CONCLUSIONS 

This paper has explored the equilibrium structure of 
an electrified interface between two phases consisting of 
charged components, as described by a phase field model. 
Such a model, being a continuum description, adds only 
the bare essentials of the physics and chemistry of electro- 
chemical interfaces: mass and volume constraints, Pois- 
son's equation, ideal solution thermodynamics in the 
bulk, and a simple description of the competing ener- 
gies in the interface. Despite this simple description, the 
model realizes the oft described behavior of the double 
layer; namely, the charge separation at the interface and 
its dependence on voltage drop (Galvani or inner poten- 
tial) across the interface. As the Galvani potential is 
varied at constant compositions of the electrode and elec- 
trolyte (constant chemical potentials), the model repro- 
duces the well known maximum of the surface free energy 
curve at the point of zero surface charge (PZC). High 
precision pseudospectral solutions of the governing equa- 
tions also deliver differential capacitance variations with 
Galvani potential that exhibit much more complex and 
realistic behavior than do the simple Gouy-Chapman- 
Stern models. The full range of behavior encompassed 
by the model must await further research. For example, 
the effect of unequal and/or nonzero barrier heights Wj 
for the components will surely affect the adsorption and, 
in turn, the surface free energy and capacitance curves. 

A recent lattice-gas model of an electrochemical sys- 
tem [24, 25] exhibits intcrfacial structures very similar 
to those found in this paper. That model also demon- 
strates simple dendrites during plating, but those lattice- 
gas papers do not explore the electrocapillary behavior 
discussed in this paper. The similarities of the predic- 
tions between that discrete model and our continuum 
approach may permit a bridge between atomistic treat- 
ments of the electrochemical interface and macroscopic 
descriptions of electroplating. 

To model a real electrochemical system with this 
method, one needs to match the parameters of the phase 
field model to the experimentally determined (or the nor- 
mally applied) understanding of the particular electro- 
chemical system. In addition to kinetic parameters de- 
scribed in Ref. [11], equilibrium solutions require several 
pieces of information. At a minimum one requires: 

1. a description of the bulk thermodynamics of the 
electrode and electrolyte, 

2. the dielectric constant of the electrolyte and elec- 
trode, 

3. an estimate for the physical thickness of the elec- 
trode / electrolyte phase interface, 

4. the actual (Galvani) potential across the interface 
for some concentration of co-existing electrolyte 
and electrode phases, and 
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FIG. 6: Comparison of the differential capacitance results of (a) this model with the predictions of (b) the Gouy-Chapman and 
(c) Gouy-Chapman-Stern sharp interface theories outlined in Appendix D and (d) the experimental measurements of Ag(100) 
electrodes in aqueous solutions of NaF [Reprinted from J. Electroanal. Chem. 138, G. Valette, "Double layer on silver single 
crystal electrodes in contact with electrolytes having anions which are slightly specifically adsorbed Part II. The (100) face", 
37-54, Copyright (1982), with permission from Elsevier]. 



5. the surface free energy or the capacitance of the 
interface for these concentrations. 

Although we currently lack an analytical expression for 
the relation between the phase field parameters K£ and 
Wj and the information in 3 and 5 above, the numerical 
results of the paper show that they are connected. In 
the future, an asymptotic analysis of the governing equa- 
tions may reveal these relationships directly. Finally, the 
non-ideal solution behavior of the electrolyte, which may 
involve complexing of ions, should be addressed. Concen- 
tration would be replaced by presumably known activity 
coefficients. 

Solution of the governing equations has proved diffi- 
cult. The resolution of the charge through the interfacial 
region requires many more mesh points than typical of 
phase field models of solidification of binary alloys. This 



is due to the more intricate structure of the charge dis- 
tribution in the interface as compared to the structure of 
the phase and concentration fields. An adaptive solution 
method, which concentrates mesh points in this inter- 
facial region, permits significantly improved calculation 
speed. 
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APPENDIX A: FIRST INTEGRAL 

Here we characterize a first integral of the steady-state 
one-dimensional equilibrium equations, and use it to ob- 
tain an expression for the surface free energy 7. 

It is convenient to introduce the electrochemical free 
energy density 



fv (£, C 1 ...C n ,<i>) = Y J Cjpj (£, Ci . . . C n , 4>) . (Al) 
j=i 

Here p,j = fXj + ZjTcj) = dfv/dCj is the electrochemical 
potential of species j, and the charge density is given by 
p = dfv/d<j). The steady-state one-dimensional equilib- 
rium equations can then be written compactly in terms 
of fv , and assume the form 



_ _Yi- _ x 



i = i 



.n- 1 
(A2a) 

(A2b) 
(A2c) 



If wc differentiate the electrochemical free energy den- 
sity of the equilibrium solution with respect to £(2), 
Cj(x), and 4>{x), we find 



- r fv(^),C 1 (x)...C n (x),cj ) (x)) 
ax 



dfvc 



(A3) 



E 

3=1 



C'Ax)-4>'(x) (e(0^) a) 



where we have used the volume constraint Eq. (4) and the 
governing equations (A2) to eliminate the partial deriva- 
tives of fv- This expression can be simplified to give 



d_ 

dx 



fv(Z(x),C 1 (x)...C n (x),<j>(x)) 
3=1 



0. (A4) 
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Since, from Eqs. (4) and (Al) we have 



f v (£(x),C 1 (x)...C n ,cf>(x)) 

n-l 

= E 



J'=l 

n-l 

J'=l 



C 3 



fJ-n 



(A5) 



[In 



we may write the first integral represented by Eq. (A4) 
in the form 



TF ~ -kH + = instant = (A6) 



<0 



V. 



v. 



where we have evaluated the integration constant in the 



far field where £ x 



and /i,, 



I'n 



In view 



of Eq. (A2a), we therefore find that the electrochemical 
potentials of the substitutional species vary through the 
interface, with 



N - l l j + V 3\-^£x 7T* 



j = l...n (16) 



The interstitial species, with Vj = 0, thus have uniform 
electrochemical potentials. 

An alternative form of the free energy functional of 
Eq. (5) takes the form 4 



F (£, Ci . . . C n 



f v ^c 1 ...c n ,4>) + ^\^- e -^-\W4>\ 2 



dV. 
(A7) 



APPENDIX B: SURFACE FREE ENERGY 

A conventional definition of the surface free energy 7 of 
a planar interface at equilibrium between two isothermal, 
multicomponcnt, fluid phases, with no electrical effects or 
volume constraints, is to write [27] 



F = J2ti™n j -P°°V + 1 A, 
3=1 



(Bl) 



where P°° is the far field value of the pressure P. The 
interface is located in the interior of the region —L/2 < 
x < L/2 and the free energy is 



F = A 



L/2 



L/2 



fv dx = A 



(B2) 



In our model of the electrolyte-electrode equilibrium, 
we are neglecting the pressure term in Eqs. (Bl) and 
(B2) , and including a volume constraint and the effects of 
an electric field on charged components. The appropriate 
definition of 7 is analogous to Eq. (Bl), with 



F 



Y^nfnj+lA., 



(B3) 



and rij is defined by Eq. (7). The surface free energy 
arises from the variation in the substitutional electro- 
chemical potentials across the interface. From Eq. (A7) 



L/2 
-L/2 



fv- 2^ c j + -rrtx - -^~4> x 



dx 



(B4) 

On substitution of Eqs. (Al) and (16) into Eq. (B4) we 
obtain 



L/2 



L/2 



" dx. 



(21) 



APPENDIX C: SURFACE CHARGE AND 
CAPACITANCE 

Here we derive the expression (19) for the variation 
in surface free energy 7 associated with changes in the 
Galvani potential A<j>° under the assumption of ideal so- 
lution thermodynamics (Eq. (25)). The derivation can 
be generalized to non-ideal solution behavior if activ- 
ity coefficients arc introduced. The variation is com- 
puted for fixed values of the far-field mole fractions X? 

and Xj, so that from Eqs. (27) and (29) we see that 
the variation SAcj) then induces corresponding varia- 
tions SA/ij = —ZjT8A(f)° and the related expression 
SfL°° = —ZjFSAcj) . The variation in A</>° also leads to 
variations 5<j)(x), SCj(x), and <5£(x) in the equilibrium 
profiles of the field variables as well. We compute the 
resulting variation of the surface free energy, £7. 
From Eq. (B4) we have 



J7 



L/2 
-L/2 



Sfv ~ E VT sc i ~ E r + K ^<^< 
3=1 3=1 



dx. 



(CI) 



4 Here we have used the identity f pc/>dx = J e(^)<p^,dx, which 
follows from the Poisson equation with appropriate boundary 
conditions. 



In computing the variation Sfv, we must consider not 
only the explicit variations arising from 6<f>(x), 8Cj(x), 
and S£(x), but also take into account the implicit varia- 
tion associated with the dependence of fv on A/i°. We 
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then find 



APPENDIX D: GOUY-CHAPMAN-STERN 



Sfv 



dfv 

d4> 



(C2) 



dfv 

on 



j=l 

0/v 



a:)-p(0p(a:)^ 



(C3) 



where we have used dfv/dAp° = Cjp(£) and SA/j,? = 
—ZjTbA<\>°. Inserting this expression into Eq. (CI) and 
integrating by parts, we find 



L/2 
-L/2 



dfv 



K^x 



e'(0 



Ete-Arw-Ew. 

J'=l 3=1 

0/v 



90 



(C4) 



Using the equilibrium equations (A2), Eq. (16), Sp,°° 
ZjT8A(f>° , and X)j=i Vj^Cj = (which follows from 
Eq. (4)), to simplify the results, we obtain 



<5 7 = <5A0° [ L/2 [l-p(£)]p(x)dx, 

J -L/2 

If we define the surface charge of the electrode as 

-L/2 

P{€)pdx 

-L/2 

and the surface charge of the electrolyte as 

-L/2 



(C5) 



(23) 



-L/2 



[1-piQpdx, 



(C6) 



Eq. (C5) recovers the classical electrochemical adsorption 
formula of Eq. (19). Note that because the total charge 



is zero, (7 



It is useful to perform a detailed comparison to the 
standard Gouy-Chapman model of the double layer. This 
model only treats variations in the electrolyte and the 
electrode-electrolyte interface is considered to be sharp. 
The inputs to the model are the difference between the 
voltage of the electrolyte at the metal 0o and the voltage 
far from the interface 0oo, the dielectric constant, and 
the cation concentration of the electrolyte far from the 
interface. The Stern modification to the Gouy-Chapman 
model requires an additional parameter x 2 -, the location 
of the plane of closest approach to the electrode of ions 
with a finite radius. The model assumes a Boltzmann 
distribution in the electrolyte and requires that Poisson's 
equation be satisfied, giving the voltage as a function of 
distance from the metal into the electrolyte, 

tanh [z M .f (</>-0oo)/4Rr] _ . _ w.gci 
tanh [z M f (fa - ^) /ART] ~ ^ [ {X Xl)l ^ J ' 

< x 2 < x < oo (Dl) 

(j) is linear for < x < x 2 and fa is the potential at x 2 
obtained by requiring continuity of cj> and of V<p at x 2 . 
The Dcbyc length of the system is 



rGC 

°<t> 

From Gauss' law 



eRT 



O/^foo 2 -7T2 



1/2 



a = — o 



dx 



(D2) 



(D3) 



and Eq. (Dl), the surface charge in the metal as a func- 
tion of <j>2 is 



too PT nl/2 • v. ZM^(02-<Aoo) 



a a = (8eC^ RT) L/ z sinh ■ 



2RT 



(D4) 



and from Eq. (20) , the differential capacitance as a func- 
tion of fa is 



1 



r> 2 T-2,/-foo \ 1/2 

RT 



cosh 



2RT 



+ *2 

e 

(D5) 



